/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "searchableSurfaceCollection.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
#include "foamTime.H"
#include "ListOps.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{

defineTypeNameAndDebug(searchableSurfaceCollection, 0);
addToRunTimeSelectionTable
(
	searchableSurface,
	searchableSurfaceCollection,
	dict
);

}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::searchableSurfaceCollection::findNearest
(
	const pointField& samples,
	scalarField& minDistSqr,
	List<pointIndexHit>& nearestInfo,
	labelList& nearestSurf
) const
{
	// Initialise
	nearestInfo.setSize(samples.size());
	nearestInfo = pointIndexHit();
	nearestSurf.setSize(samples.size());
	nearestSurf = -1;

	List<pointIndexHit> hitInfo(samples.size());

	const scalarField localMinDistSqr(samples.size(), GREAT);

	forAll(subGeom_, surfI)
	{
		subGeom_[surfI].findNearest
		(
			cmptDivide  // Transform then divide
			(
				transform_[surfI].localPosition(samples),
				scale_[surfI]
			),
			localMinDistSqr,
			hitInfo
		);

		forAll(hitInfo, pointI)
		{
			if (hitInfo[pointI].hit())
			{
				// Rework back into global coordinate sys. Multiply then
				// transform
				point globalPt = transform_[surfI].globalPosition
				(
					cmptMultiply
					(
					    hitInfo[pointI].rawPoint(),
					    scale_[surfI]
					)
				);

				scalar distSqr = magSqr(globalPt - samples[pointI]);

				if (distSqr < minDistSqr[pointI])
				{
					minDistSqr[pointI] = distSqr;
					nearestInfo[pointI].setPoint(globalPt);
					nearestInfo[pointI].setHit();
					nearestInfo[pointI].setIndex
					(
					    hitInfo[pointI].index()
					  + indexOffset_[surfI]
					);
					nearestSurf[pointI] = surfI;
				}
			}
		}
	}
}


// Sort hits into per-surface bins. Misses are rejected. Maintains map back
// to position
void Foam::searchableSurfaceCollection::sortHits
(
	const List<pointIndexHit>& info,
	List<List<pointIndexHit> >& surfInfo,
	labelListList& infoMap
) const
{
	// Count hits per surface.
	labelList nHits(subGeom_.size(), 0);

	forAll(info, pointI)
	{
		if (info[pointI].hit())
		{
			label index = info[pointI].index();
			label surfI = findLower(indexOffset_, index+1);
			nHits[surfI]++;
		}
	}

	// Per surface the hit
	surfInfo.setSize(subGeom_.size());
	// Per surface the original position
	infoMap.setSize(subGeom_.size());

	forAll(surfInfo, surfI)
	{
		surfInfo[surfI].setSize(nHits[surfI]);
		infoMap[surfI].setSize(nHits[surfI]);
	}
	nHits = 0;

	forAll(info, pointI)
	{
		if (info[pointI].hit())
		{
			label index = info[pointI].index();
			label surfI = findLower(indexOffset_, index+1);

			// Store for correct surface and adapt indices back to local
			// ones
			label localI = nHits[surfI]++;
			surfInfo[surfI][localI] = pointIndexHit
			(
				info[pointI].hit(),
				info[pointI].rawPoint(),
				index-indexOffset_[surfI]
			);
			infoMap[surfI][localI] = pointI;
		}
	}
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::searchableSurfaceCollection::searchableSurfaceCollection
(
	const IOobject& io,
	const dictionary& dict
)
:
	searchableSurface(io),
	instance_(dict.size()),
	scale_(dict.size()),
	transform_(dict.size()),
	subGeom_(dict.size()),
	mergeSubRegions_(dict.lookup("mergeSubRegions")),
	indexOffset_(dict.size()+1)
{
	Info<< "SearchableCollection : " << name() << endl;

	label surfI = 0;
	label startIndex = 0;
	forAllConstIter(dictionary, dict, iter)
	{
		if (dict.isDict(iter().keyword()))
		{
			instance_[surfI] = iter().keyword();

			const dictionary& subDict = dict.subDict(instance_[surfI]);

			scale_[surfI] = subDict.lookup("scale");
			transform_.set
			(
				surfI,
				coordinateSystem::New
				(
					"",
					subDict.subDict("transform")
				)
			);

			const word subGeomName(subDict.lookup("surface"));
			//Pout<< "Trying to find " << subGeomName << endl;

			const searchableSurface& s =
				io.db().lookupObject<searchableSurface>(subGeomName);

			// I don't know yet how to handle the globalSize combined with
			// regionOffset. Would cause non-consecutive indices locally
			// if all indices offset by globalSize() of the local region...
			if (s.size() != s.globalSize())
			{
				FatalErrorIn
				(
					"searchableSurfaceCollection::searchableSurfaceCollection"
					"(const IOobject&, const dictionary&)"
				)   << "Cannot use a distributed surface in a collection."
					<< exit(FatalError);
			}

			subGeom_.set(surfI, &const_cast<searchableSurface&>(s));

			indexOffset_[surfI] = startIndex;
			startIndex += subGeom_[surfI].size();

			Info<< "    instance : " << instance_[surfI] << endl;
			Info<< "    surface  : " << s.name() << endl;
			Info<< "    scale    : " << scale_[surfI] << endl;
			Info<< "    coordsys : " << transform_[surfI] << endl;

			surfI++;
		}
	}
	indexOffset_[surfI] = startIndex;

	instance_.setSize(surfI);
	scale_.setSize(surfI);
	transform_.setSize(surfI);
	subGeom_.setSize(surfI);
	indexOffset_.setSize(surfI+1);
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::searchableSurfaceCollection::~searchableSurfaceCollection()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
{
	if (regions_.size() == 0)
	{
		regionOffset_.setSize(subGeom_.size());

		DynamicList<word> allRegions;
		forAll(subGeom_, surfI)
		{
			regionOffset_[surfI] = allRegions.size();

			if (mergeSubRegions_)
			{
				// Single name regardless how many regions subsurface has
				allRegions.append(instance_[surfI] + "_" + Foam::name(surfI));
			}
			else
			{
				const wordList& subRegions = subGeom_[surfI].regions();

				forAll(subRegions, i)
				{
					allRegions.append(instance_[surfI] + "_" + subRegions[i]);
				}
			}
		}
		regions_.transfer(allRegions.shrink());
	}
	return regions_;
}


Foam::label Foam::searchableSurfaceCollection::size() const
{
	return indexOffset_[indexOffset_.size()-1];
}


Foam::pointField Foam::searchableSurfaceCollection::coordinates() const
{
	// Get overall size
	pointField coords(size());

	// Append individual coordinates
	label coordI = 0;

	forAll(subGeom_, surfI)
	{
		const pointField subCoords = subGeom_[surfI].coordinates();

		forAll(subCoords, i)
		{
			coords[coordI++] = transform_[surfI].globalPosition
			(
				cmptMultiply
				(
					subCoords[i],
					scale_[surfI]
				)
			);
		}
	}

	return coords;
}


void Foam::searchableSurfaceCollection::findNearest
(
	const pointField& samples,
	const scalarField& nearestDistSqr,
	List<pointIndexHit>& nearestInfo
) const
{
	// How to scale distance?
	scalarField minDistSqr(nearestDistSqr);

	labelList nearestSurf;
	findNearest
	(
		samples,
		minDistSqr,
		nearestInfo,
		nearestSurf
	);
}


void Foam::searchableSurfaceCollection::findLine
(
	const pointField& start,
	const pointField& end,
	List<pointIndexHit>& info
) const
{
	info.setSize(start.size());
	info = pointIndexHit();

	// Current nearest (to start) intersection
	pointField nearest(end);

	List<pointIndexHit> hitInfo(start.size());

	forAll(subGeom_, surfI)
	{
		// Starting point
		tmp<pointField> e0 = cmptDivide
		(
			transform_[surfI].localPosition
			(
				start
			),
			scale_[surfI]
		);

		// Current best end point
		tmp<pointField> e1 = cmptDivide
		(
			transform_[surfI].localPosition
			(
				nearest
			),
			scale_[surfI]
		);

		subGeom_[surfI].findLine(e0, e1, hitInfo);

		forAll(hitInfo, pointI)
		{
			if (hitInfo[pointI].hit())
			{
				// Transform back to global coordinate sys.
				nearest[pointI] = transform_[surfI].globalPosition
				(
					cmptMultiply
					(
					    hitInfo[pointI].rawPoint(),
					    scale_[surfI]
					)
				);
				info[pointI] = hitInfo[pointI];
				info[pointI].rawPoint() = nearest[pointI];
				info[pointI].setIndex
				(
					hitInfo[pointI].index()
				  + indexOffset_[surfI]
				);
			}
		}
	}


	// Debug check
	if (false)
	{
		forAll(info, pointI)
		{
			if (info[pointI].hit())
			{
				vector n(end[pointI] - start[pointI]);
				scalar magN = mag(n);

				if (magN > SMALL)
				{
					n /= mag(n);

					scalar s = ((info[pointI].rawPoint()-start[pointI])&n);

					if (s < 0 || s > 1)
					{
					    FatalErrorIn
					    (
					        "searchableSurfaceCollection::findLine(..)"
					    )   << "point:" << info[pointI]
					        << " s:" << s
					        << " outside vector "
					        << " start:" << start[pointI]
					        << " end:" << end[pointI]
					        << abort(FatalError);
					}
				}
			}
		}
	}
}


void Foam::searchableSurfaceCollection::findLineAny
(
	const pointField& start,
	const pointField& end,
	List<pointIndexHit>& info
) const
{
	// To be done ...
	findLine(start, end, info);
}


void Foam::searchableSurfaceCollection::findLineAll
(
	const pointField& start,
	const pointField& end,
	List<List<pointIndexHit> >& info
) const
{
	// To be done. Assume for now only one intersection.
	List<pointIndexHit> nearestInfo;
	findLine(start, end, nearestInfo);

	info.setSize(start.size());
	forAll(info, pointI)
	{
		if (nearestInfo[pointI].hit())
		{
			info[pointI].setSize(1);
			info[pointI][0] = nearestInfo[pointI];
		}
		else
		{
			info[pointI].clear();
		}
	}
}


void Foam::searchableSurfaceCollection::getRegion
(
	const List<pointIndexHit>& info,
	labelList& region
) const
{
	if (subGeom_.size() == 0)
	{}
	else if (subGeom_.size() == 1)
	{
		if (mergeSubRegions_)
		{
			region.setSize(info.size());
			region = regionOffset_[0];
		}
		else
		{
			subGeom_[0].getRegion(info, region);
		}
	}
	else
	{
		// Multiple surfaces. Sort by surface.

		// Per surface the hit
		List<List<pointIndexHit> > surfInfo;
		// Per surface the original position
		List<List<label> > infoMap;
		sortHits(info, surfInfo, infoMap);

		region.setSize(info.size());
		region = -1;

		// Do region tests

		if (mergeSubRegions_)
		{
			// Actually no need for surfInfo. Just take region for surface.
			forAll(infoMap, surfI)
			{
				const labelList& map = infoMap[surfI];
				forAll(map, i)
				{
					region[map[i]] = regionOffset_[surfI];
				}
			}
		}
		else
		{
			forAll(infoMap, surfI)
			{
				labelList surfRegion;
				subGeom_[surfI].getRegion(surfInfo[surfI], surfRegion);

				const labelList& map = infoMap[surfI];
				forAll(map, i)
				{
					region[map[i]] = regionOffset_[surfI] + surfRegion[i];
				}
			}
		}
	}
}


void Foam::searchableSurfaceCollection::getNormal
(
	const List<pointIndexHit>& info,
	vectorField& normal
) const
{
	if (subGeom_.size() == 0)
	{}
	else if (subGeom_.size() == 1)
	{
		subGeom_[0].getNormal(info, normal);
	}
	else
	{
		// Multiple surfaces. Sort by surface.

		// Per surface the hit
		List<List<pointIndexHit> > surfInfo;
		// Per surface the original position
		List<List<label> > infoMap;
		sortHits(info, surfInfo, infoMap);

		normal.setSize(info.size());

		// Do region tests
		forAll(surfInfo, surfI)
		{
			vectorField surfNormal;
			subGeom_[surfI].getNormal(surfInfo[surfI], surfNormal);

			const labelList& map = infoMap[surfI];
			forAll(map, i)
			{
				normal[map[i]] = surfNormal[i];
			}
		}
	}
}


void Foam::searchableSurfaceCollection::getVolumeType
(
	const pointField& points,
	List<volumeType>& volType
) const
{
	FatalErrorIn
	(
		"searchableSurfaceCollection::getVolumeType(const pointField&"
		", List<volumeType>&) const"
	)   << "Volume type not supported for collection."
		<< exit(FatalError);
}


void Foam::searchableSurfaceCollection::distribute
(
	const List<treeBoundBox>& bbs,
	const bool keepNonLocal,
	autoPtr<mapDistribute>& faceMap,
	autoPtr<mapDistribute>& pointMap
)
{
	forAll(subGeom_, surfI)
	{
		// Note:Tranform the bounding boxes? Something like
		// pointField bbPoints =
		// cmptDivide
		// (
		//     transform_[surfI].localPosition
		//     (
		//         bbs[i].points()
		//     ),
		//     scale_[surfI]
		// );
		// treeBoundBox newBb(bbPoints);

		// Note: what to do with faceMap, pointMap from multiple surfaces?
		subGeom_[surfI].distribute
		(
			bbs,
			keepNonLocal,
			faceMap,
			pointMap
		);
	}
}


void Foam::searchableSurfaceCollection::setField(const labelList& values)
{
	forAll(subGeom_, surfI)
	{
		subGeom_[surfI].setField
		(
			static_cast<const labelList&>
			(
				SubList<label>
				(
					values,
					subGeom_[surfI].size(),
					indexOffset_[surfI]
				)
			)
		);
	}
}


void Foam::searchableSurfaceCollection::getField
(
	const List<pointIndexHit>& info,
	labelList& values
) const
{
	if (subGeom_.size() == 0)
	{}
	else if (subGeom_.size() == 1)
	{
		subGeom_[0].getField(info, values);
	}
	else
	{
		// Multiple surfaces. Sort by surface.

		// Per surface the hit
		List<List<pointIndexHit> > surfInfo;
		// Per surface the original position
		List<List<label> > infoMap;
		sortHits(info, surfInfo, infoMap);

		// Do surface tests
		forAll(surfInfo, surfI)
		{
			labelList surfValues;
			subGeom_[surfI].getField(surfInfo[surfI], surfValues);

			if (surfValues.size())
			{
				// Size values only when we have a surface that supports it.
				values.setSize(info.size());

				const labelList& map = infoMap[surfI];
				forAll(map, i)
				{
					values[map[i]] = surfValues[i];
				}
			}
		}
	}
}


// ************************************************************************* //
